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ABSTRACT 

A system of two partial differential equations represents the 
transient heat transfer behavior of compact heat exchanger surfaces 
when subjected to a step change in fluid temperature. A solution is 
presented for this system which includes the effects of longitudinal 
thermal heat conduction. Also presented are the solutions for the two 
limiting cases of zero and infinite longitudinal conduction. The 
numerical results have been compared to those of C. P. Howard indicating 
a significant decrease in computational time and an increase in accuracy 
of results, The revisea curves of maximum slope of fluid temperature 
versus NTU should be of practical value in the evaluation of heat- 
transfer data obtained by transient testing of compact heat exchanger 
surfaces. An unusual combination of mathematical techniques is presented 
for the solution of a boundary value problem involving partial differen- 
tial equations. The solution combines the application of Laplace trans- 
formation with a numerical technique eu eoned by Ho Hurwitz in. and 
P. F. Sweifel, and adapted by L. A. Schmittroth for the inversion of 
Laplace transforms. This technique greatly expands the number of cases 


to which Laplace transforms may be successfully applied. 
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NOMENCLATURE 
Z 

total heat transfer area, ft 

; 7 
matrix cross sectional area, ft . 
step increase in fluid temperature, = Vie oe 
specific heat of fluid, BTU/1b °F. 
specific heat of solid, BTU/1b °F. 


2 on 


heat transfer coefficient, BTU/hr ft 
thermal conductivity of matrix, BTU/hr ft °F. 


total length of fluid flow path, ft. 


dimensionless conduction parameter, = KA,/wc 


dimensionless heat transfer units, = hA/w ce. 
dimensionless time variable, = hAO [WI co. 
time, hr. 

sole temperature, °F. 

reference solid temperature, °F. 

fluid temperature, °F. 

fluid temperature at X = 0, t = 0 +; °F. 
reference fluid temperature, °F. 

fluid mass flow rate, lb/hr. 

mass of matrix, 1b. 


distance along the fluid flow path, ft. 


dimensionless position variable, =X 2. 


f° 





1. Introduction. 

The “single blow" problem refers to the study of the transient heat 
transfer behavior of compact heat exchanger surfaces for the purpose of 
determining their heat transfer characteristics. Interest in compact 
heat exchangers has grown with the development of gas turbine engines. 

A regenerative cycle is desirable to increase the efficiency of the gas 
turbine. Proper design of a compact heat exchanger requires knowledge 
of the heat transfer characteristics of the various materials in numer- 
ous geometric configurations. Until recently it has been necessary to 
design and build the complete heat exchanger in order to determine its 
heat transfer characteristics experimentally. The objectives of the 
single blow studies are to determine the heat transfer characteristics 
from small test sections of the various matrices and to obtain consider- 
able reduction in experimental costs. 

The aim of the experimental testing is to obtain the exit fluid time- 
temperature history subsequent to a step change in the entrance fluid 
temperature. The experimental history must be compared to a theoretical 
time-temperature history to determine NTU, the dimensionless heat trans- 
fer parameter. The heat transfer coefficient, h, or the Colburn J factor 
can then be calculated. 

In 1950 Locke [9] outlined five experimental eecnniauee). One of the 
simpler of these has become known as the'"maximum slope" method. This 
technique compares the maximum slope of the experimental time-temperature 
history to the theoretical maximum slopes for various values of NTU and 
of sv , the longitudinal conduction parameter. With a cursory inspection, 
one might not appreciate the efficiency of this single point curve matching 


ayambers in square brackets refer to the bibliography. 
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technique. One maximum slope represents a complete theoretical or a 
complete experimental curve. Therefore, it is unnecessary to calculate 
and plot vast quantities of theoretical time-temperature histories, 

This method eliminates the laborious and inaccurate matching of theoreti- 
cal and experimental time-temperature histories, previously required to 
determine NTU. However, the maximum slope technique has two major limita- 
tions. The first problem is the difficulty in analytically solving the 
system of governing differential equations when longitudinal conductivi- 
ties other than zero or infinity are included. The second problem is the 
magnification »f experimental error due to curve matching. This problem 
is discussed in detail later. 

Considerable analytical effort has been made to solve various aspects 
of the single blow problem since 1927 when Nusselt first studied it. 
Schumann [14 | developed a solution in 1929 for zero conductivity in the 
direction of flow. His results are for a porous medium such as gravel but 
have been extended to include matrices consisting of metal balls, wire 
screens, or continuous materials for which the effects of longitudinal 
thermal conduction can be ignored. However, as indicated by Mondt [11 | and 
Creswick [3 | , usually when the matrix is constructed of a continuous 
material in the flow direction, the effects of longitudinal conduction 
must be considered. Creswick outlined a finite difference technique to 
include this effect but his work was not complete enough to cover the area 
encountered in experimental testing. Mondt EY obtained a closed solu- 
tion for the limiting case of 7A = © . Figure 1 indicates how drasti- 


cally these two cases differ and shows the need for intermediate curves. 
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Maximum Slope 
Figure l 

In 1948, F. E. Romie, et. al. [12 | developed a system of partial 
differential equations governing the transient heat transfer behavior for 
hollow circular cylinders. Creswick arrived at the same set of equations 
for parallel plates. By a finite difference technique, he numerically 
approximated the solution of the system. However, difficulty with con- 
vergence was sometimes encountered. In 1963, C. P. Howard [5 | » guided 
by physical considerations, developed an alternate form of the finite- 
difference equations, and he was able to determine the convergence criteria 
and avoid the problem encountered by Creswick. Howard presented a table 
of results and a set of curves of NIU versus Maximum Slope for various 
values of A . The table covers a range quite adequate for most experi- 
mental testing. As the values of NTU increased beyond 60, it became very 
difficult to obtain solutions with the finite difference technique due to 
the large computation times involved. 

The objective of the solution presented in this paper is to verify or 
improve the results obtained by C. P. Howard and to avoid large computation 


times as NTU increases. This solution not only avoids the difficulty at 





large values of NTU, but becomes faster as NTU increases due to the 
corresponding increase in the time at which the maximum slope occurs. 
Referring to the error magnification due to curve matching, C. P. 
Howard [5] showed in Figure 3-A of his Appendix 3, that the error due 
to the experimental determination of maximum slope is magnified by a 
factor of two or greater when the value of NIU is determined by the 
maximum slope technique. The error magnification factor is defined as 
the ratio of percent error in NTU to percent error in experimental 
Maximum Slope. The essential features of the error factor as a function 


of NTU are indicated in Figure 2, 
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This figure shows that for the special case, fl = 0, the lowest 
error magnification is between 2.5 and 2 for NTU greater than five. At 
NTU near two there is a peak for which an accurate determination of NTU 


by maximum slope matching is impossible. When A is non-zero,the peak 





shifts to slightly higher values of NTU and the useful range above NTU 
of four decreases as f\ increases. For example, for A = .06, an 
error factor of 3, at NTU = 5, is doubled at NTU = 28. Therefore, with 
this technique, one is forced to accept much larger errors as NTU and 

fA increase. To avoid the inherent magnification of error, a modifica- 
tion can be made to the time-temperature curve maiching technique with 
the aid of modern computers. When q_ is finite, non-zero it is manda- 
tory to use a computer to calculate the theoretical values of fluid 
temperature, given NTU, and C It would be elementary to include a 
least Square error curve matching technique in the temperature program. 
Then, given three or more experimental values of fluid temperature and 
their corresponding values of t, the computer would compare the calcu- 
lated theoretical temperatures at the given t's for various NTU's and 
determine the NTU with the least square error. The anticipated error 


in any particular theoretical value of temperature is + 1% or less. 


z : 
"t'™ denotes the dimensionless time parameter. 





tags Mathematical Technique. 

The approach to be presented solved the system of partial differ- 
ential equations for finite, non-zero A _ by eliminating the variable 
of solid temperature. The result is a third order partial differen- 
tial equation of fluid temperature as a function of time and x (the 
ratio of distance along the matrix to the total length). This equa- 
tion and the necessary boundary conditions were transformed with re- 
spect to time by Laplace transformations to obtain a third order total 
differential equation with respect to x with parameter s. The equation 
was then solved for the transformed fluid temperature. The transformed 
boundary conditions were applied to determine the three arbitrary co- 
efficients which were functions of the parameter s. Because the trans- 
formed solution was very complicated, numerical inversion was used to 
compute the inverse Laplace transform, 

A Gaussian quadrature method, developed by H. Hurwitz, Jr. and P. F. 
Sweifel [6] for Fourier Transform Integrals, was adapted by L. A. Schmitt- 
roth [13] for the numerical inversion of Laplace Transforms. This techni- 
que avoids the difficulty encountered with alternating, slowly converg- 
ing functions and proved to be essential for the successful inversion of 
this solution. For the two limiting cases A = 0 and A= , a direct 


inversion of the Laplace solution was available. 





3. Solution, 
The governing differential equations for the transient heat transfer 


behavior including a finite, non-negative, longitudinal conductivity are 


as potions 


(1) OMe = (nr—u) + ea 
NTU OX 


aU 
Sy OU NTU (u-rv) 


The applicable boundary conditions are the following: 








a. V(x,0) = 0 d. UU (0,t) =0 
Ox 

b, OM (x,0) = 0 e. seltk (1,t) = 0 
Ox OX 


€: ( VQOste =.6 


The following is the approach used to solve this set of equations 


OW 
ot 


respectively. In a similar manner the solid temperature, u, and other 


bf 


for the fluid temperature, v, and the slope of fluid temperature, 


desired values can be obtained. 


Solving equation 2 for u results in 


(2 22 See eelde 
ge NTU ax OT 


The partial derivative of u with respect to t and the second partial of 








u with respect to x are determined from 2a, and then substituted into 


equation 1 to give 








3 2 
_ oft) qe OA-(x,t) _ NTU OAT (x,t) 
Zz 
__ NTU ONT (x,t) NTU OM (X,t) 
nr ot Xot 


>the equations were developed from a heat balance by Creswick[3] ; 





Each term in equation 3 is transformed with respect to t by Laplace 
transformation along with the following boundary conditions: 


a. v(x,0) = 0 








els 
Be ,0) = 
7X (x, 0) 
The necessary transforms a 
Aah) = 2 es 
Cc, re ae = Bx. ( S) soe. 
al ise 
d. Lj 2 Gt) =): rete = aoe | +S [ardxt Je dt 
= —ar(%0) + SLar(x3) | 
Baru, of — Sar 5 _ 
= a I CX,S ) 
ae ea 0) +5 


Substitution of these values in equation 3 after applying the above bound- 


ary conditions, results in 


_ _ 9(XS) NTU a os) NTU(S +1) O(*.S)_ NTU S V5) 


The corresponding auxiliary equation is 


Z 
(3b). me + NTU AE — STO C544) ee in =o 


The general solution in the Laplace s-plane for the fluid temperature is 


ny X N2X IL 3X 


4). ar(y,s)= Cie +C62s)€ + C3(S)€ 


where rl, r2, r3 are the roots of equation (3b). 
The boundary conditions c, d, and e are transformed and then used to 


determine the coefficients Cl, C2, and C3. The transform of C is 


al? 
3) -— —— 
UV koe ) 





The transform of d is QUES) = OQ , and similarly for e, 


2 
Cree CX, tae ole 5 (X,¢) 
O X 
Applying boundary condition ¢ to equation 4 gives 


(6). PAO) = Ome rates) + C3Co) = G 
Ss 
Rather than apply boundary conditions d and e directly it is convenient to 
use their equivalent form in terms of v. To apply boundary condition d 


it is necessary to transform equation 5 which results in 


2 
au (x,s)_ 4_ 29M (x,5)_ gar(y,S) therefore; 
ax NTU OX2 By 


4 
CGS) eee 1UL CO) 2) an (0,5) — 
aX — NTU OX* ae Px O 


When this boundary condition is applied the following equation is obtained: 


(7). A (recy Sry Ce 4575 C5(8)) + nO, +*n2CeS)+N3 C; GCG) S07: 


The application of boundary condition e is similar, resulting in 


ny h 
(8). A (ntce ens Ce ee ise Janie reese Nz (,e°°= 0, 


NTU 
Zz 
(8a). Define Rm= on + Tm) miecemneet 2, 3, Then ene 


following matrix equation can be set up using equations 6, 7, 8. 





G 

Ss 4 Ne un GAS) 
0 = Ri Re Rs ; C2 CS) 
O Re” Re Re’ | 3) 


To solve for Cl, C2, C3, Cramer's Rule is applied. The denominator is 


A= [RR:(eve™ )+ R Ro(e™-e™ ) + RiR2(e™—E™ )| , 


The three arbitrary coefficients are the following: 


CZ 1 4 
O Ro R3 
X n 
Cisy—- 10 Ree* Rew] _ ela Re Rare” | 
A A : 
‘ee 
Ry O Ra and 
e | 
C,(S) = Re" O Poe. me eS RRC" _RPe™ | 
A IS. 
tL A Os 
nur Ro O 
cya tRe™ Rem o | _SURRae™ Ree" 
IS IS 


Substituting Cl, C2, and C3 into equation 4 results in 
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© _ ¢ [RRet-e™ Je™+R Re" Je™ RR(eN e" Ee 
uOs) = = nN 


Evaluating equation 9 at x = 1 results in the following: 


(40) 
r). ath (RitMa) — Chs+N2) + (7,+113) 
eo cS Reper he Japp frre Jn arm e] 
5 ReRs(e" _e” )+R(e™-e" ) + RRe(e™ —e” ) 


Recall that the Laplace operation Sf(Ss) — F (+0) corresponds to 
d F(t). It follows directly that av (Ae) is equal to 
dt at 

the operation Se LS). , since ach 0) = GED, 


boundary condition QA . Therefore, 


QI) Bar (4,5) = 


or Sear (4,8) 


To determine the exit fluid temperature and the slope with respect 
to time of the exit fluid temperature, it is necessary to find the in- 
verse Laplace transforms of equations 9 and 10 respectively. If the 
values of R, and res i=1, 2, 3; were not functions of s this step would 
be elementary. Unfortunately, the values of R and r are indeed functions 
of s as indicated by equations 8a and 3b respectively. Therefore, for 
this particular case the task of finding an analytic inverse transform 
is nearly hopeless. 

Rather than proceed further on this approach, a computer program to 
calculate the complex roots for given values of s, NTU, and 7A , is 
used. Consequently, this step precludes any possibility of finding an 


analytical solution if one exists. Now, the value of v(1,s) or our at 


at 


(1,s) can be determined for any set of the parameters, C, s, NIU, and ri : 


ll 





and the Laplace inversion integral may be used to find v(l,t). Find- 
ing an adequate numerical technique to evaluate this inversion integral 
was no simple matter. 

When Simpson's Rule was applied to evaluate the inversion integral 
using the real (VR) and imaginary (VI) parts of v(l,s), the results were 
grossly in error, the behavior being as described in reference |6] and [13}. 
The numerical technique described in reference [13] was found to be greatly 
superior to Simpson's Rule for this particular problem. It incorporates 
a technique to accelerate convergence of an alternating series, and uses 
Gauss-Chebyshev formulas for integration. This technique uses either 
. the real part (VR) or the imagirary part (VI) and results in two differ- 
ent forms. Since there are no known values for this solution it was 
necessary to use both of these forms to be confident that the correct 
solutions were obtained. 


According to the theory of Laplace inversion, the inversion integral, 


Ei J. top YU £%S) eC ds ,» should be independent of G, where 


U 
s =G+ iy . When this approach was used with Simpson's Rule none of 
the results were independent of G. Recall the form of the inversion 
integral = f WG, 9) otyt dy . Since the region 
initially investigated used t = 7 and G of .5, 1 and 2, any error in the 
numerical technique was greatly multiplied giving very poor results, However; 
this indicated that G is an additional parameter that must be determined 
for each new program. Rather than cause difficulty, this parameter proves 
to be essential. For the same G, it was found that one program gave con- 


sistently low results while the other was consistently high. Then by 


“(Elncsebvisl Op A"elte ep. ery OF 
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varying G of both programs the program using VR can be made to agree 
with the results of the program using VI. This is the approach that 
was used to make the various programs converge to a solution. The 
numerical technique outlined in [13] uses Chebyshev's polynomial to fit 
the functior f(s). It was found that this polynomial is very sensi- 
tive to a small change in G. For example, upon changing G from 0.0 to 
0.05, the product (slope e NTU) changed from 1.0877 to 1.1333 an increase 
of 4.19%. Again according to theory, results should be independent of G. 
Therefore, two new programs were written using a Gauss-Legendre quadra- 
ture format to integrate VR or VI. These programs were found to be 
essentially independent of G for reasonable increments of G. For example 
changing G from 0.0 to 0.06, the slopeeNTU changed from 1.10096 to 1.09776 
a decrease of 0.29%. Figure 3 indicates a comparison of these two pro- 
grams. On this basis alone, it is felt that the Legendre polynomial is 
more accurate. Unfortunately, in comparison to the programs using the 
Chebyshev polynomial it is considerably slower. Because of this fact, 
the Legendre programs which include an error analysis are being used to 
determine accurate test values for the determination of the best G for 
each program. Once the best G is determined, the speedier Chebyshev 
programs can be used for long data runs. 

The integration technique described in detail in references [6 | and 
[13] will be described here as it was applied to this problem. The inver- 


sion integral is the following: 


3 Gteoo 
“(re =f) = 4 { £6) e “ds 
G 





21 L Dives 


Then writing f(s) as the sum of the real and imaginary parts results in 
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FIGURE 3 








i ff [ P(oriy) Hoty) a! dy) 


The inversion oe can then be written as 


(12a) Cip=£ jf Love ign ieoriselty + {fs [perignin(orihe9" dy) 


where the following condition applies: f(s) zs ay. 





It follows directly that 
(12b) DC6+Ly) =P9(G-iy) , and 
(12e) 44 (G+iy) = — 4(@-iy), 


In the first integral of equation 12a, let y =-y} and dy = -dy', then 


CO —Ly't : 
(13) £()=2< 1 f [D(G-iy') + We~ivledy' 4 if Fp(Gsiy) si4l(oxcy)|e sy) 


Applying the relations 12b and 12c, equation 13 reduces to 


flys = [4 ‘Ld oriy)-i4icoriy)| et “+[douy)+iacoriy)]e dy 


On rearranging and simplifying, this becomes 


Ge aD ze 
(14) £(4) = ef | }Co+iy) cos Wil HG +iy) sinyt | dy 


By definition of the inverse integral, f(-t) = 0. Thus 


O= a" | $(ercy) COS yt +ucoriy) singt] dy 


EE 


5 
A bar over a quantity indicates the complex conjugate. 
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which reduces to 
OD B 0O . : t 
J dcor+iy) cos ytdy = ~f yl(eriy) Singtady. 
O 
Two forms of f(t) can be found by expressing the integrals either in 


terms of the real or the imaginary parts. They are respectively 


Gt ,CO 
(16) F(t)= SS { Dloriy)cosyt dy and 


Gt 
Qn) Li) BE yi(riy) Sin yt dy. 


These two equations are the basis of the two programs mentioned earlier. 
In general discussions VR will represent the function P cos yt and VI 
the function 4! sin yt. The Gaussian quadrature formula using either 
Chebyshev's or Legendre's polynomial is used to evaluate the integral 
from zero to infinity. Both functions VR and VI are of the general 


form shown in Figure 4 


or 
Wale 


Figure 4 
If the function were integrated with respect to y and the summation 


stopped at a, there would be significant error. Two operations are 
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included in this solution to reduce the error. The first operation 
is an averaging technique which accelerates convergence. The second 
operation involves integrating each half wave and then summing them. 
To accomplish this we introduce the following change of variable in 


Equation 17. 


Let * = = y- £ therefore y= (x+ 5) and dy= Tox, 


Applying these to equation 17 results in 


(i7a) fF(t)= — ee Naan (X+4)]Sintt(x + Yr) OX , 
= aes GO | 
eee =. =f [+i xt E)] sintr(x+h)dy, 


Now let x X—- Mm , then 


Gt 


Q7) £(ty= — = 5 ee UT (xX +eme4)|Sintt (x+m+2qx’ 


M=0 


Upon supressing the prime, we have 


(18) £(¢)— +iZ (xtm+)|COS 1X Ox 


~ 
o, 
Au 
re \ 
= 

oe 


When equation 18 is written in terms of partial sums, the result is 


ay Fe) = 2S Ty 





ta mn—» 9 where 


Den (t)= _ alps (t) and 





Ge i= ey U1] G ri (emedy) COS (1TX) Ax : 


Again referring to tas L. A. Schmittroth applied a Gaussian quadrature 


method to evaluate T(t). IT, xt) can b approximated as follows: 


4 oN N 
(19b) i Goes (-4)" “> Wi (Hwy + Hees | 
jet SEIS I 
whe re 
(19c) yw = 2i-4 We ONO, RN 


Z2(2N+4) 


N 
a  G) — EC Y 40+ 4) 5 Pee (4; +m+s), 


The number of points used to fit the function is 2 N, Obviously, the 
larger N is, the greater the degree of accuracy one can expect. This 
is the core of the iterative process and an adequate N is desired to en- 
sure reasonable computer time. This point will be discussed in detail 


later. The N coefficients Ww are the solutions of the system. 


(19d) 7? x W; Nagg*h | @j-4) 7 | 2 A L +z) 
J=A 2(2N +4) {1 ' | '(e+4) 
whe re b sie hea ee a 
Equation 19 uses Chebyshev's polynomial for the curve fitting. A 
similar form of io t) can be developed from Equation 16. As mentioned 
earlier, an integration technique was developed using a Legendre poly- 
nomial for curve fitting. If equation 16 is used to illustrate this 


approach, the ralation a = ty is required. Then equation 16 reduces to: 
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Ge prs | 
(20) f(t) pie’ D(6+i EL) Cos 7X dx. 


Now, the Gaussian quadrature method can be applied as above. For the 
Gauss-Legendre formulas, a change of variable was applied to equation 17b 
to set up the infinite sum of integrals from zero to one, one to two and 
so on rather than a sum of integrals evaluated from -1/2 to 1/2. 

The Gauss-Chebyshev and the Gauss-Legendre formulas applied to 
the real or imaginary form of the inversion integral produce four differ- 
ent computer programs. Any one of these should give the correct solution, 
but one form may prove to be better for a particular application. 

For zero and infinite longitudnal conduction, the governing differ- 
ential equations reduce to much simpler forms. The solutions for both 
cases are given in Appendix I and are in complete agreement with the 


solutions of Schumman and Mondt for A = 0 and A = OO respectively. 
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4, Computer Program, 

The principal objective of the computer program is to solve for the 
maximum ree, fluid temperature at X = Le, The main program is design- 
ed to determine the slope of the fluid temperature at X =fwhich is the 
numerical inversion of equation ll. 

The subroutine SOOTS2 calculates the real and imaginary coefficients 
of the equation 3b given the parameters NTU, A, G, and Y. These co- 
efficients are then put into subroutine TOOTS2. 

Subroutine TOOTS2 calculates the real and imaginary parts of the 


roots r r, of equation 3b. It sends its results back to SOOTS2 


ie 2a 
which arranges them in the proper order and subscripts them appropriately. 
SOOTS2 then provides the proper roots as input values to EVAL. The EVAL 
Picnutiine provides the main program with the values of the real (VR) 
and imaginary part (VI) of a i S> » given the parameters G, 
TTX, C, and the roots from SOOTS2. 

The main program will call for N of these values (see equation 19c) 
with which it will calculate T(t) (see equation 19). The programmer is 

wn 

free to determine the number of I (t)'s to sum, where Z I m(t)= 555 (t+) - 
Fifteen to twenty partial sums appear to be sufficient for the function en- 
countered in this problem. Each partial sum calculated is stored and 
serves aS an input to subroutine AVER. 

The purpose of subroutine AVER is to accelerate convergence of the 

partial sums by an averaging technique. This technique is discussed in 
reference | 6| - The subroutine will take an arbitrary number of partial 


e h ° 
sums and calculate any specified n° average, and return it to the main 


th 
program. It also calculates the difference between successive n averages 
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returning it to the main program, The main program then requires this 
difference to be less than an arbitrarily specifiec value (EP). If the 
difference is larger than EP, the main program increases the number of 
partial sums by LD and calls AVER again until the difference is less 
than EP. At this point, the value of the accelerated partial sum is 
accepted and the slope is calculated. However, the slope calculated for 
a particular value of dimensionless time is not the value desired. It 
is necessary to search for the maximum slope by varying time. 

Initially, a search using an incremental step of time was used. 
This was eventually abandoned because it would find relative maxima only 
and it was tediously slow if a poor choice were made for a starting time. 
It should be noted that theoretical relative maxima have been found for 
small values of NTU and negative slopes have also been noted. J. M. 
Bannon[1] is concurrently conducting an experimental thesis and has experi- 
mentally found relative maxima and negative slopes for high mass flow or 
small values of NTU. It would be interesting to compare a theoretical 
time temperature curve for his experimental value of A with his ex- 
perimental time-temperature curve. 

The purpose of SLOMAX is to find the maximum slope. To avoid the 
problem of choosing a poor starting time, a reasonable range of t is 
chosen, A rule of thumb is that the largest maximum slope will occur at 
a time less than or equal to NTU. SLOMAX uses a parabola curve fitting 
technique with three points. 

The first and last of the three points are program input values of 
time and the middle point is the arithmetic mean. The slopes are calcu- 


lated for the first two values of time and then tested to ensure that the 
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middle slope is larger than the first. This is necessary to obtain the 
maximum rather than a minimum slope. Then the third value of slope is 
calculated and reared to insure that it is less than the middle slope. 

The program has the ability to adjust the arbitrary values of time to 
ensure the middle slope is the largest. Once this basic shape is obtained 
the program calculates the value of time at which the parabola has a maximum. 
This process is repeated using the calculated maximum time and slope as the 
middle point and the two nearest points from the preceeding iteration. 

When the difference of the new value of maximum slope and the previous 
calculated maximum slope is less than EPl, a program input, the slope is 
accepted as the maximum. Then the values of time, slope and input para- 


meters are printed out. 
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TABLE 5. ‘RESULTS OF SLO2 AND SLO4 FOR = 0.00 
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Or. Discussion of Results. 

Table 1 indicates that programs SLO3 and SLO5 are not sensitive to 
G and that both programs give the same solution to six decimal places 
for NTU greater than 5. The variation of G from -0.1 to +0.1 is much 
larger than the variations used in the SLO2 and SLO4 programs which are 
more sensitive to G. All results presented were calculated with eight 
point formulas. 

One way to increase the accuracy of the program is to increase the 
number of terms in the partial sum. There is a test in the program 
which will increment the number of partial sums if the nth average has 
not converged to within 1 x ine of the preceeding nth average for each 
calculation of slope. The test for the sufficient number of partial 
sums indicates that for NTU greater than five, twenty partial sums will 
give the same accuracy as 100 partial sums for programs SLO3 and SLOS5. 
At NTU of five, the SLO3 program automatically increased the number of 
partial sums to 25 which then gave as accurate a solution as 95 partial 
sums. There is an indication that for lower values of NIU, more partial 
sums may be necessary. 

Rew. Maxim [10] presented a table of values of maximum slopes cal- 
culated with Schumann's solution for the special case of A = 0.0. 
Schumann's solution is an infinite series of Bessel functions. Each 
value of slope was calculated using Bessel functions accurate to ten 
significant places and the summation was continued until there was no 
change in the tenth decimal place of the solution. Unfortunately, the 
search by Maxim could not distinguish between relative maximum slope 


values and his search ended when the slope at the next increment of time 
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was less than the preceeding slope. Relative maximum slopes have been 
found by the theoretical solution presented in this paper and their ex- 
istence is supported by experimental evidence found by J. M. Bannon i 1] ; 
Figure 26. The results of Table 2 clearly indicate that Maxim's values, 
though accurate to six places, were calculated at a time much too low to 
have found the largest relative maximum slope in the NTU range less than 
five. For values of NTU equal to and greater than five, both the time 

at which maximum slope occurred and the values of maximum slope times NTU 
calculated by SLO3 agree to at least seven decimal places with Maxim's 
results. SLO5 solutions were in complete agreement with the results of 
SLO3. 

The computation time depends on the degree of accuracy required to 
calculate one value of slope. The accuracy required of the search routine 
for the results of Table 2 was + 1 x 10°” difference in the new value of 
maximum slope as compared to the previous iteration. For values of NTU 
less than een, the number of iterations required was four or five times 
the number for NIU values greater than ten. Once the approximate value 
of time at maximum slope is known, the search routine can be set to find 
the maximum with fewer iterations. But, even without optimum starting 
times the computation time was approximately half a minute for the A= 0.0 
case, For non zero values of A , the computation time per search is 
expected to be twelve minutes, 

Results of SLO5S for yr. =.04 are shown to be insensitive to G in 
Table 3. One test point at NIU = 10 indicates good agreement between 
SLOS5 and SLO3 programs. As for the jt = 0 case, results indicate good 
agreement with C. P. Howard's results [2] » For NIU greater than five 
Howard's results are consistantly low but all have less than 1% differ- 


ence as compared to SLO5. As in the A = 0 case, Howard's results for 
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NTU equal to two are much lower than those obtained by SLO5. Again this 
is probably a failure of his search technique to distinguish between 
relative maximum slopes. Unfortunately, he did not publish the times 

at which the maximum slopes occurred. Howard also used the incremental 
step search and the search ended when the newly calculated slope was 
less than the preceeding slope. 

Programs SLO2 and SLO4 are quite sensitive to parameter G. The 
results obtained graphically and listed in Table 4 were compared to the 
accurate solution of Maxim for AU = 0.0. The percent error of the graphi- 
cal results was less than + 1.5% except at NTU = 2. Figure 6 indicates 
that this intercept was not determined. Considerable effort would be 
required to find this intercept and obviously the mean value at NTU = 2 
has unacceptably large error. For the other values of NTU, however, 
intercepts were easily found as indicated by Figure 5. 

Table 5 shows similar results for A = 0.005. No accurately known 
values are available for this region. C. P. Howard's results are com- 
pared to the results obtained by SLO2 and SLO4. If as expected, the 
SLO2 and SLO4 results are 1 to 1.5% too large, then C. P. Howard's 
solutions for JA = 0.005 would be approximately .5% too low for NTU 
greater than 10. Figure 7 again shows how the intercept or best G was 


obtained. 
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car Conclusions and Recommendations. 

Program SLO5, using the VR, Gauss-Legendre formulas, is a fast 
accurate tool capable of supplying all the theoretical maximum slope 
data needed for experimental testing of compact heat exchanger surfaces, 
The SLO3, VI, program is equally accurate but slightly slower. The SLOS 
program is being applied to verify the NTU versus Maximum Slope curves 
presented by C. P. Howard [5] » The results of these calculations will 
be published at a later date. 

The combination of Laplace transforms with numerical inversion should 
be applicable to most boundary value problems in which one is able to 
express the solution as an inverse Laplace transform. Therefore, this 
technique is primarily limited in application to linear systems of part- 
ial differential equations with constant coefficients. 

The two main problems encountered in the numerical inversion were 
the slow convergence of the integration and radical behavior of the func- 
tion of s at time equal zero. The VR and VI functions were cyclic in 
nature with slowly decreasing magnitudes and period. For this particular 
solution the Legendre polynomial approximated the functions better than 
Chebyshev's polynomial. However, for other functions a number of other 
polynomials may give better results. The acceleration of convergence is 
adequately accomplished by the present numerical inversion but it gener- 
ates a problem in addition. Because of the averaging technique applied 
to the partial sums to accelerate convergence the formal error analysis 
normally applied to each partial sum is impossible. Therefore, the only 
way remaining to determine the accuracy of the calculated maximum slopes 


was by comparison with R. E. Maxim's accurately determined values for the 
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limiting case of fA = 0. Even though direct inversion of the Laplace 
transform was possible for this case, the numerical inversion was ac- 
complished in seconds with the same accuracy attained by the approxi- 
mation to an infinite sum of Bessel functions carried out to six decimal 
place accuracy. For non-zero jf , there were no accurate solutions with 
which to compare. For these values an indication of accuracy was obtain- 
ed by calculating the same points with an independent program. The SLO3 
program is not independent of SLO5 in a strict interpretation of indepen- 
dence. However, they are different programs, SLO3 using the imaginary 
part of the inversion integral and SLO5 using the real part. Either 
program should give the same result. The same result to six or seven 
decimal places was obtained with these programs indicating excellent 
accuracy for a numerical process. The comparison of SLO5 results with 
Howard's results, obtained with a rather crude finite-difference tech- 
nique, indicated that he obtained very accurate results by judicious ap- 
plication of the finite-difference technique. Most of his results were 
less than one half of one percent low when compared to SLOS. 

The following recommendations indicate avenues of investigation 
that could be undertaken with either computer programs developed in con- 
junction with this thesis or the mathematical technique presented. 

It is recommended that the temperature programs TEMP4 and TEMP2 be 
applied to obtain a time-temperature history for comparison with corres- 
ponding experimental results to determine how well the mathematical model 
represents the behavior of the experimental test rig. 

An investigation should be conducted to determine the cause of and 
physical significance of the relative maximum slopes obtained both experi- 


mentally and theoretically for the heating cycle. One objective would be 


Sp; 


to determine if the relative maximum obtained experimentally correspond 
in magnitude and time of occurance with those obtained by theory, Ex- 
perimentally it was noted that the relative maximums are practically 
indistinguishable for the cooling of the solid, but that for the same 
matrix and mass flow of fluid the relative maximums are very distinct 
for the case of heating the solid. An investigation should be made to 
determine theoretically the effect of longitudinal conduction for the 
cooling cycle as compared to the present solution for the heating cycle. 
The objective of this investigation would be to determine if the theoreti- 
cal solution is capable of predicting a difference in response on cool- 
ing as compared to heating. The cooling cycle solution would require 
appropriate sign changes in the heat balance from which governing partial 
differential equations were derived. This may, therefore, require a new 
solution for cooling since the sign changes could cause a change in the 
auxiliary cubic equation 3b of Section 3. It is possible that a Sign 
change of a term in this equation could cause a change in the significance 
of the parameter A on the result of fluid temperature or slope of fluid 
temperature. Bannon [i | has noted experimentally that the largest maximum 
slope obtained by cooling agreed fairly well with that obtained by heating. 
The theoretical results, in the region of NTU = 2, are difficult to 
obtain, and the determination of NTU given an experimental maximum slope 
in this NTU region is impossible. An attempt should be made to develop 
a curve matching technique employing a least square error method. This 
method would use three or more values of time and their corresponding 
temperatures from an experimental time-temperature history. The TEMP4 


program would calculate theoretical values of temperature at the given 
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times for an estimated value of NTU. Then an NIU search routine could 
calculate the sum of the squared error of each temperature. After re- 
peating this operation for two bracketing values of NIU, a Lagrange 
iteration method would calculate the NTU with the least squared error. 
Another possible means of curve matching would be to match the dimension- 
less time at which the maximum slope occurs for the region of NTU = 2, 
rather than matching the magnitudes of maximum slope. This would require 
an adequate starting mark on the experimental time-temperature history 
and a correction of the time at the maximum slope to account for the 
system time delay of the experimental response. 

A more general system of governing partial differential equations 
can be obtained from the heat balance as derived bs Creswick, by in- 
cluding a term to represent the energy stored in the fluid. This term 
is required for transient heat transfer testing using a liquid rather 
than a gas. The result of this addition changes equation 2 of section 


3 to the following 


OW A, [u-w- ae S| 
y} 


OX NTU ot 
where n¢ we eCeAg ek , the dimensionless fluid capacitance, 
Ws Cs 


and Y is the density of the fluid in Tyas, This equation combined 
with Equation 1 of Section 3 would then be the new governing partial 
differential equations. The present mathematical technique should be 


capable of solving this system of equations. 
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APPENDIX I 


Solution for Zero and Infinite Longitudinal Conduction. 


1. ZERO LONGITUDINAL CONDUCTIVITY 


As in the solution for finite conductivity the development of this 


case is the samc initally; cf, Equation 3a, page 8. 


(1) dr (xs) partes). (s+) u(x 5) _ US ar(%s) = =0 
3x3 Ox 
it is necessary to multiply by a. 


For this case where A is defined as a, 


The substitution of a = 0 results in the equation 


(2) ae SS. Wwe= © 
oe 


The general solution of equation 2 is 


—-bS x 
G@ nri,s)= Cte" = Cie 


The boundary conditions are the same as in the previous solution 


Boundary conditions <j and b were applied in order to arrive at equation 


C which transforms into v(0,s) = 


1, above. Boundary condition C is v(0,t) = 


C/s. Applying aay condition ¢ to equation 3 results in 





nr (0,S)= Cc - C\(s) , therefore 
oS 

4 (V = ~ 6S x nd i C are 
" (X,5)= le ee 3 vhs) == Ee 


S 
It follows directly, 
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AW (x oe 3 (1,5) polos 
(5) Ot Ce aX for ot = (1c Sri 


This special case was solved analytically in 1929 by T. E. W. 
Schumann [14] . The following development will show that the present 


solution agrees with Schumann's. 





Let s' =s + 1 and apply it to equation 4. 
/ 
=f 
Then -b2 : X or 
mix, s-t={7C E 5 
ee 


©) ar(x,S-t) = Cen*¥s et 


COs 


4 
sim 


From page 229 of reference [2 | , we find 


M PCA e&) _ ft = 
Suakes S — j;_—~ 
sm ey Ludi (2. VKt ) 
Comparing equation 6 and 7, the inverse transform of equation 6 becomes: 


oo es ‘farwst} = Ce "s (E = Th. (2 VER). 


From page 294 [2] is 

-L at ~4 
i { £(s-a)§ = © Pe SFiS) therefore, 
(9) W(X,t)= Cefe™ ‘= (= tye * T,(2-V6ex) 


(9a) CX,t) = ce bOE) F ( ()2 ia (21eex ) 
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Define a = btx and W = 2 -fa’ , then compare equation 9a with the form 


4 
(9b) a2 T,(2va) | 
Note that 
{VY} 
<\* aa m off 
(9c) (S —— = ie 6. 
ox (bxt) “2 


Equation 9a then becomes 


~ beat Oma 


M=0 


(10) W(x t) = Ce 


Schumann's solution is 


(11) ar(y,zZ) 2C °° Ss sal Ala Je (2i-Syz) 


M20 d(yz)” : 


= Tn (2Ve), 


Equations 10 and 11 are equivalent which can be shown as follows. From 


page 392 [4] , the modified Bessel function identify is 


fw" 1, (WN) 2) TE ew) 
ad Ww 
Then, for n = 0 and W = 2 7/a, 


—— ie = 

o(W) a (W 
Ty; | 4 ) 
Apply the chain rule such that 


GEG aly 

da IW da 
It follows that | 
GIT, ave = @* JT (2va) 
da | 


where M=O 





By mathematical induction it can be shown that 
MY 
mM payee 
(12) boxe =Q 2 [m(z{e), 
= 


Equation 10 then becomes 


= i 
(13) wx, t) = Ce is a td [I.(zva)| 
m=O ao 
From page 392 [4] ’ 


I m(21@) = ("In(20 VO"), 
T, (ata) = Jy (2cve"). 


Note that Schumann's development is for the case of heating the fluid as 


and 


time increases and that the present case is developed for cooling of the 
fluid with increasing time, which changes the sign of t in the exponential 


term . It then follows that equation 13 is, indeed, 


a) mriyz)= ee Fy zd” [To(eeVue)] 
(ee Ce 
Since the program for numerical inversion of equations 5 and 6 is already 
available it is much easier and probably much faster to apply it than to 
evaluate equation 11. The same theoretical solution was obtained by G. L. 


Locke in 1950 [9] : 
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s* ' 2 = ———> => —=_ 
= > one @&—» 
== 


ities ee we MV =| ee i — it ; 


eee st sett) om eeecee cee) eed aad 


Bit Mire [2b Py E 32.87 JME 


Tat 







eveezpie ot) bee f aw le mee ‘elec leer se) eee ale eee 
b) GA) OF viret  'serey tc PMS Bo Okey ee ay ef OORT 
4 oo UP hetthey, Geer ney ener ao; 64) ep anepe 


th orst & 


2. INFINITE LONGITUDINAL CONDUCTIVITY 


When A = © , the solid temperature is no longer a function of x. 


a“ u (x,t) = 
In fact the following is true: yKm ) =o Therefore, equation 


1, page 7, becomes indeterminate and must be replaced by:/ 
1 
OU 
(1) Se Ct) i [ Ar (x,t) -U(t)| Ax and equation 2 is 
o 


oe -~ NTU | ut) -ar(x,t)] 


The boundary conditions to be applied to this case are: 
a. v(x,0) = 0, 


b. v(x,0) 0, and 


ey ey. (06) @= 967) 6 
Solve equation 2 for u(t) and take the derivative of u with respect 
to t. Then substituting these back into equation 1 and applying boundary 


condition c, results in 


= 
1 aut), Bar (4%) = [ fore = )—nr(x,t)| dx 


NTU oxot 
= = at all 


Multiply through by NIU and rearrange the terms to arrive at 


a IATOG) 4 Tu _dar (xt) Parlin = C= © 
oxot dt 


Transformation of each term in 3 by Laplace gives 


or further discussion, see page 47. 
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(4) Jar (x0), Sdar(x,S) a nTU[-mcx, 0) +Sarcx,s)] +nr(t,5) —S 20 
uae may +N J Ss 


Applying boundary conditions a and b and dividing by S$, equation 4 reduces 


to 


(5) Or (K,S) NS iO ae C _ WAS) , 
dX 55 S 
Again by treating s as a parameter,equation 5 is a total differential 


equation. The solution of the differential equation is 


(6) ~ArlxX.s) = Cr(sy OCT = .£ [C _ arts) 
( 9 ) Gr Se ata Va & S ( 


Now apply boundary condition c to equation 6 which results in 





NTU S | 


(7) RnCOnS) . = Gs) + 4 Sie . - (4s) 





Therefore, Cy Cue aE _ AL | 4 a (Ape SY. and 
< NTu -S NTU +S , 





(8) _ 0 SCE) C _ po NTU+X 
Vixns)= =e 2 va as 3] lS ) 


Solve equation 8 for v(l,s* as follows: 


pL SY) = @ oN +t.[¢- v(h3)] (\-e""™ ) 


ees ; pace J] Z ‘a cane 


NTu - S © NTU - S 
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- NTU —~NTU 
Define ( (4 - = ) ; then €@ =e" yee , and 


The particular fluid temperature of interest is that temperature at x = l, 
Also the value of C must be dimensionless, for comparison with Mondt's 


solution, and can be set equal to one. It follows directly that 


(MUP pS ures ae — $-NTU) 


(S+ ¢) 


a aras)2 —tL ~ £eTY , Sf 
Sia S+¢ Sct 


The inverse Laplace transform of equation 10 is 


~ Se -Sé ae 
Wy (4,t) =C — NTU-SE +f{-e which reduces to 
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oe 
Q1) ardét,t) = 1- NTU Se 


y) 
~NTU 
le 
where & = ap » The solution as a function of 
time is Sar 
Geer) 





Rec) = ie = eC NTU t 

This solution is in complete agreement with that derived by J. R. 
Mondt [11] in 1961. Since a convenient inverse transform was available 
it was unnecessary to apply the intricate numerical inversion. 

To derive Equation (1) for this special case it was necessary to 
refer to the heat balance from which Creswick [3 | derived the governing 
equations. The last term of Creswick's Equation (1) was eliminated and 
the remaining terms were arranged in terms of dimensionless parameters. 
Then both sides of the equation were integrated over the length of the 


matrix to arrive at Equation (1) on page 44. 
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APPENDIX II 
Flow diagrams and program listings. 
The flow diagrams will preceed the applicable program listings. 


The symbols used in the flow diagrams are defined as follows: 


Terminal - the beginning, end, or point of 
interuption in the program 


Input/Output - input symbols are preceeded by 
a terminal. 


Decision - branches are labeled by sign 
according to the sign of the 
value enclosed, 

Connector - enclosed numbers refer to the 
prefix numbers of the corres- 
ponding statements on the pro- 
gram listing. 


< | Offpage Connector - designates an entry to or exit 
from a page, 


There are six distinct programs which are used separately, although 
they could be combined into a larger package to make use of one or another 
at the user's option. There are the four ''slope'’ programs, SLO2, SLO3, SLO4 
and SLO5 which have been described in detail in the body of the thesis, and 
the two "temperature" programs TEMP2 and TEMP4, which have been referred to 
in the text but not described in detail. The latter are intended to pro- 
vide.a time-temperature history; no use of these programs has been made in 


this thesis. See page 35 for recommendations concerning their use. 
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PROGRAM SLO3* 





Set 11=0 











Cee ely 
Ose LO! SS 
3 L1, LD, It data card 














+ 
= 
\ (Output) (Output) 
Evaluation Uses Legendre 
subroutine Polynomial 
2 Short 













Do 203 J=1.,100. 
Call GaussQ | 


2 
= 


a 
eke 
21) 2 I1, Slope 


Call Szomax | 
Tii4=8IM/3 


*Also applies to program SLO5 









_ L=L+LD (IX 


‘Call Ave ey ~ 
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PROGRAM SLo2” 


: \G, B, A, T1,-TM Set I1 =0 
(107) TD, EP, L1 LD, Cc) On last 


11 data card 







= <a> | 
“Pp 
Commu) (Output) 
Uses Evaluation 
Chebyshev Subroutine 
Polynomia 25 
a 











200 
500 SU(I)= 
0.0 





{Call SLOMAX 
TM1 = 8T™ 





DO 203 J=1, 100 


603 DO 200 I= bo 8 
Y = PI*p 
CALL SOOTS2: 


CALL EVAL2S | “Gor 
*Also applies to SLO4 ~ Ho : = : 
2000 | | 


= P1+(1./9.) 















COuupur) 
SUM, DIF 
SuU(J) 


(Output) 
I1, Slope 
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PROGRAM TEMP" 


6) BP e}-O- SE 
™, t 


D, EB, On last - 
Ieee D. ola data card 
M1 
SEop 10 (10) ; <> 
— 


















(Cucpu ce 
Program is 
Tempera - 


(Output) 
Evaluation 
Subroutine 


1 Short olynomial 





Do 203 ciate 100 
Call cs 


L=L?TLD 
Caiat 
AVER 
Gc} > (203 
gp A 


fe) 
| Continue 


T=T+TD =a 
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